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ABSTRACT 

T-H ! 

O 

^^ ' We present a series of simple, largely analytical models to compute the effects of 

disruption on the mass function of star clusters. Our calculations include evaporation 
^ • by two-body relaxation and gravitational shocks and mass loss by stellar evolution. We 

find that, for a wide variety of initial conditions, the mass function develops a turnover 
"^ . or peak and that, after 12 Gyr, this is remarkably close to the observed peak for globular 

^T) . clusters, at Mp ss 2 x 10^ Mq. Below the peak, the evolution is dominated by two- 

^ I body relaxation, and the mass function always develops a tail of the form ip{M) = 

Q"^ I const, reflecting that the masses of tidally limited clusters decrease linearly with time 

^j I just before they are destroyed. This also agrees well with the empirical mass function 

(•~^ ' of globular clusters in the Milky Way. Above the peak, the evolution is dominated 

by stellar evolution at early times and by gravitational shocks at late times. These 

processes shift the mass function to lower masses while nearly preserving its shape. The 



o 

:^. 

O , . radial variation of the mass function within a galaxy depends on the initial position- 

velocity distribution of the clusters. We find that some radial anisotropy in the initial 
H I velocity distribution, especially when this increases outward, is needed to account for 

52 I the observed near-uniformity of the mass functions of globular clusters. This may be 

consistent with the observed near-isotropy of the present velocity distributions because 
clusters on elongated orbits are preferentially destroyed. These results are based on 
models with static, spherical galactic potentials. We point out that there would be even 
more radial mixing of the orbits and hence more uniformity of the mass function if the 
galactic potentials were time-dependent and/or non-spherical. 

Subject headings: celestial mechanics, stellar dynamics — Galaxy: kinematics and dy- 
namics — galaxies: kinematics and dynamics — galaxies: star clusters — globular 
clusters: general 
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1. Introduction 

Globular clusters appear to have a preferred mass scale. Their mass function has a turnover 
or peak at Mp ~ 2 x 10^ Mq and a dispersion of only cr(logM) ~ 0.5 (Harris 1991). The mass 
and luminosity functions of globular clusters are often modeled as lognormal functions, although 
they can also be represented by several broken power laws (McLaughlin 1994). Moreover, globular 
clusters represent the most numerous gravitationally bound stellar subsystems within the spheroidal 
components of galaxies (the others being dwarf galaxies). There are relatively few subsystems in 
galactic spheroids with masses between those of individual stars and globular clusters (i.e., in the 
range lO" ^ M < 10^ -^o) a-^d between those of globular clusters and the spheroids themselves 
(i.e., in the range 10^ ^ M < lO^'' Mq for a large galaxy). The preferred mass scale of old 
globular clusters may be an important clue in understanding their formation and evolution. The 
corresponding feature in the luminosity function, at My ~ —7.3, is sometimes used as a distance 
indicator (e.g., Whitmore 1997). 

In contrast, the mass functions of many other types of astronomical objects appear to be scale- 
free and are often modeled as single power laws. For diffuse and molecular clouds in the Milky Way, 
the mass function has the form ^(M) oc M^ with /3 ~ —2 over the range 10~^ ^ M < 10^ Mq 
(Dickey &: Garwood 1989; Solomon & Rivolo 1989). For young star clusters in the disks of normal 
galaxies (i.e., open clusters), the luminosity functions, which may or may not reflect the mass 
functions, are power laws, (t>{L) oc L" with a ~ —2 (Milky Way: van den Bergh & Lafontaine 1984; 
Large Magellanic Cloud: Elson k Fall 1985; M33: Christian & Schommer 1988). Finally, for the 
young star clusters formed in interacting and merging galaxies, the luminosity and mass functions 
also have power-law form, with a f» /3 ~ —2 for 10*^ ^ M < 10^ Mq (Whitmore et al. 1999; Zhang 
& Fall 1999, and references therein ^). The last finding is particularly significant because these 
clusters are often regarded as young globular clusters. 

Two explanations have been proposed for the preferred mass scale of old globular clusters. One 
is that the conditions in ancient galaxies and protogalaxies favored the formation of objects with 
masses ~ 10^-10^ Mq but that these conditions no longer prevail in modern galaxies. For example, 
the minimum mass of newly formed star clusters, set by the Jeans mass of interstellar clouds, will 
be high when the gas cannot cool efficiently and low when it can, which in turn will depend on the 
abundances of heavy elements and molecules, the strength of any heat sources, and so forth. These 
effects may have suppressed the formation of low-mass clusters in the past but not at present (Fall 
& Rees 1985; Kang et al. 1990). The other explanation for the preferred mass scale of old globular 



'^Whitmore et al. (1999) found that a double power law with ai = —1.7 and 02 — —2.6 provided a slightly 
better fit than a single power law with a ~ — 2 to the luminosity function of young star clusters in the Antennae 
galaxies (NGC 4038/9). However, this still contrasts sharply with the luminosity function of old globular clusters. 
Suggestions that the luminosity and mass functions of the young star clusters in the Antennae might be similar to 
those of old globular clusters were based on earlier, less sensitive observations, which could not reach beyond the 
putative turnovers in these functions (Meurer 1995; Fritze-von Alvensleben 1999). 



clusters is that they were born with a much wider spectrum of masses that was later modified 
by the selective destruction of low-mass clusters (Fall &: Rees 1977; Gnedin & Ostriker 1997, and 
references therein). In this case, a power-law mass function might evolve into a lognormal-like mass 
function (Vesperini 1997, 1998; Baumgardt 1998). This idea is appealing because the masses and 
sizes of the brightest young clusters in merging galaxies are similar to those of the old globular 
clusters in the spheroids of galaxies. 

Star clusters are relatively weakly bound objects and are vulnerable to disruption by a variety 
of processes that operate on different timescales. Stellar evolutionary processes remove mass from 
clusters by a combination of supernovae, stellar winds, and other ejecta. These are effective on 
both short timescales (t < lO'^ yr), when the clusters or protoclusters are partly gaseous, and 
on intermediate timescales (lO'^ ^ t < few x 10* yr), when the clusters consist entirely of stars. 
Three stellar dynamical processes remove mass from clusters on long timescales (t > few x 10* yr). 
First, internal relaxation by two-body scattering causes some stars to gain enough energy to escape 
from the clusters. Second, as clusters orbit around a galaxy, they experience a time-dependent tidal 
field, which may vary rapidly enough when they pass near the galactic bulge or through the galactic 
disk that stars in the outer parts of the clusters cannot respond adiabatically. The corresponding 
changes in the energy of the stars (heating and relaxation) cause some of them to escape. These 
effects are known respectively as bulge and disk shocks and more generically as gravitational shocks. 
Third, dynamical friction, the deceleration of clusters induced by the wakes of field stars or dark 
matter particles behind them, causes the clusters to spiral toward the galactic center, where they 
may be destroyed by the strong tidal field. 

A potentially serious problem with the idea that disruption causes the turnover in the mass 
function of globular clusters is that the chief disruptive processes operate at different rates in 
different parts and different types of galaxies (Caputo & Castellani 1984; Chernoff, Kochanek, & 
Shapiro 1986; Chernoff & Shapiro 1987; Aguilar, Hut, & Ostriker 1988; Gnedin & Ostriker 1997; 
Murali & Weinberg 1997a, b,c). For example, the rate at which stars escape by two-body relaxation 
depends on the density of a cluster, which is determined by the tidal field, and hence is higher in 
the inner parts of galaxies than in the outer parts. The rate at which stars escape by gravitational 
shocks is also higher in the inner parts of galaxies, both because the orbital periods are shorter there 
and because the surface density of the disk is higher there. Moreover, disks are absent in elliptical 
galaxies. Thus, if the mass function were strongly affected by disruptive processes, one might 
expect its form to depend on radius within a galaxy and to vary from one galaxy to another. This, 
however, appears to be contradicted by many observations showing that the luminosity function of 
clusters (a mirror of the mass function when the spread in ages is relatively small) varies little, if 
at all, within and among galaxies (Harris 1991). 

The goal of this paper is to explore the evolution of the mass function of star clusters by 
a variety of disruptive processes, including evaporation by two-body relaxation and gravitational 
shocks and mass loss by stellar evolution. We are especially interested in how the mass function is 
affected by different position- velocity distributions of the clusters, and which of these are compatible 



with observations. We formulate this problem in terms of a simple, approximate model that can 
be solved largely analytically. This clarifies how the mass function is affected by the different 
disruptive processes and different position-velocity distributions. Our calculations are performed 
in the context of static, spherical galactic potentials. But we also discuss qualitatively how our 
results would be affected by time-dependent and non-spherical galactic potentials. The plan for 
the remainder of the paper is as follows. In Section 2, we specify our model, with the associated 
assumptions, equations, and parameters. We present the results of our calculations in Section 3, 
showing the influence of each physical effect on the evolution of the mass function. In Section 4, 
we compare our results with previous studies, and in Section 5, we summarize our conclusions. 



2. Models 

We are interested here in the evolution of the mass function of star clusters, defined such that 
ip{M, t)dM is the number of clusters with masses between M and M + dM at time t. We assume 
that the clusters present initially, at t = 0, lose mass continuously (smoothing over the escape 
of individual stars) and that no clusters are created subsequently. Then the mass function must 
satisfy the continuity equation 

where the dot denotes differentiation with respect to t. The formal solution of this equation is 

V'(M,t) = V'o(Mo)|aMo/aM|, (2) 

where iPq{M) = ^(M, 0) is the initial mass function, and Mo(M, t) is the initial mass of a cluster 
that has a mass M at a later time t. In addition to these variables, the mass function may depend 
on the orbital parameters of the clusters, and hence their location within a galaxy, and may also 
vary from one galaxy to another. We first consider clusters on the same orbit. Later, we will 
average the mass function over realistic distributions of orbits and examine its dependence on the 
properties of the host galaxy. 

We consider three processes that reduce the masses of star clusters: evaporation driven by 
two-body relaxation, evaporation driven by gravitational shocks, and mass loss driven by stellar 
evolution (supernovae, stellar winds, and other ejecta). Dynamical friction, combined with tidal 
limitation, also causes disruption, but this is only important near the centers of galaxies and 
is neglected here, mainly to simplify our analysis. More specifically, a cluster of mass M at a 
distance R from the center of a galaxy with a circular velocity Vc would be destroyed in a time 
tdf ~ 60{Vc/220 km s"^)(M/10^ MQ)-'^{R/kpcf Gyr (see equation 7-26 of Binney & Tremaine 
1987). We have evaluated this expression for the 146 globular clusters in the Milky Way with known 
luminosities (assuming M/Ly = 3) and positions in the most recent compilation of data by Harris 
(1996, 1999). We find that only two clusters have tdf < 10 Gyr and only seven have tdf < 20 Gyr; 
the vast majority have much larger t^f and are thus virtually immune to disruption by dynamical 



friction. This suggests, but does not prove, that dynamical friction was also relatively unimportant 
in the past. Bulge shocks contribute to the disruption of clusters on highly elongated orbits, but 
they are probably less important than disk shocks for most clusters and are also neglected here. 
The rates of evaporation by bulge and disk shocks depend on the properties of the clusters and 
their orbits in similar ways. Thus, by including strong disk shocks near the centers of galaxies, we 
also mimic at least qualitatively the effects of bulge shocks. 

As an approximation, we assume that the processes considered here — two-body relaxation, 
gravitational shocks, and stellar evolution — operate independently of each other, at fractional rates 
Uev, ^sh-, and Use- Thus, following many previous studies (see Spitzer 1987), we write 

M = -{yev + ysh + Vse)M, (3) 

with 

trh M'^/'^rl' 
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In equation (4), ^g is the fraction of stars that escape per half- mass relaxation time t^h by two-body 
scattering, r^ is the half-mass radius, m is the mean stellar mass, and In A is the Coulomb logarithm. 
We neglect possible slow variations in the last two quantities and set m = 0.7 Mq and InA = 12 
in all our calculations. In equation (5), tgh is the gravitational shock heating time for first-order 
energy changes in the impulse approximation (Ostriker, Spitzer, & Chevalier 1972) and ^ is a 
correction for partial adiabatic (i.e., non-impulsive) response averaged over all the stars in a cluster 
(see Appendix A for details). The factor 7/3 accounts approximately for the addition of second- 
order energy changes, also known as shock-induced relaxation (Spitzer & Chevalier 1973; Kundic & 
Ostriker 1995). The other coefficient in equation (5) relates the fractional change in energy caused 
by gravitational shocks to the corresponding fractional change in mass, i.e., M jM = HgE/ E. Also 
in equation (5), Vz is the vertical component of the velocity of a cluster relative to the galactic 
disk, P^ is the azimuthal period of its orbit around the galaxy, and gm = 2ttGTi(1 is the maximum 
vertical acceleration caused by the disk of surface mass density S^. We assume that the disk has an 
exponential profile all the way into the galactic center, T,(i{R) = T,(i{0)exp(—R/R(i), thus helping 
to mimic the effects of bulge shocks. The fractional rate of mass loss by stellar evolution depends 
on the age of a cluster and the stellar initial mass function (IMF). We compute I'seit) from the 
Leitherer et al. (1999) models with the Salpeter IMF. 

We assume each cluster has an outer, limiting radius rt determined by the tidal field of the host 
galaxy at the pericenter of its orbit. Clusters on orbits with fixed pericenters, such as those in a 
static, spherical galactic potential, as assumed here, will therefore evolve at constant mean density, 
p = M/(47rrf/3); and clusters on orbits with different pericenters will have different p (with an 
additional weak dependence on the shape of the orbits, resulting from the centrifugal acceleration 
at pericenter; see equation (15) below). This is a standard assumption, although it is not expected 



to be perfect except possibly for circular orbits (Spitzer 1987). The assumption that the clusters 
are tidally limited is justified by the fact that, if they initially extended beyond rt, their outer 
parts would be stripped off after a few orbits, whereas if they did not initially extend to rt, they 
would expand as a result of the disruptive effects considered here, predominantly stellar mass loss 
in the early stages, until they reached r^. Some clusters with low central concentrations might be 
destroyed relatively quickly by a combination of stellar evolution and tidal limitation, with little or 
no help from two-body relaxation or gravitational shocks (Chernoff & Weinberg 1990; Fukushige 
& Heggie 1995); these clusters are not included in our calculations. 

We must now specify the escape probability parameter for two-body relaxation ^g, the energy- 
mass conversion factor for gravitational shocks Kg, and the relation between the half- mass and 
tidal radii, r/j and rt- A valuable point of reference is Henon's (1961) model for the self-similar 
evolution of a tidally limited cluster with a single stellar mass by two-body relaxation alone. This 
has ^e = 0.045 and rh = 0.145ri. The Henon model is often regarded as an adequate approximation 
for high-concentration clusters before core collapse, which typically occurs about half way through 
their lifetimes, and an excellent approximation for all clusters after core collapse. The value of ^e 
found in Monte Carlo and Fokker Planck models with a single stellar mass is typically 2-3 times 
below the Henon value in the early stages of evolution and closer to it in the late stages (Spitzer & 
Chevalier 1973; Lee & Ostriker 1987; Gnedin, Lee, & Ostriker 1999). On the other hand, models 
with a realistic spectrum of stellar masses evolve a few times faster than those with a single stellar 
mass (Johnstone 1993; Lee &: Goodman 1995). Thus, in most of our calculations, we adopt the 
Henon value of ^e as a reasonable approximation to the effective escape probability parameter for 
the entire evolution of a realistic cluster, including both its pre- and post-core collapse phases. We 
also adopt the relation between rh and rt in the Henon model. 

The energy-mass conversion factor Kg depends on how the energy imparted to a cluster by grav- 
itational shocks is divided between bound and escaping stars. The detailed, but non-evolutionary 
calculations by Chernoff et al. (1986) give k^ ~ 1 for high-concentration clusters (k^ is —vjf in 
their notation). In the evolving Monte Carlo models of Spitzer &; Chevalier (1973), which include 
two-body relaxation and the first- and second-order energy changes caused by impulsive gravita- 
tional shocks, the total evaporation rate is given to an accuracy of about 30% by equation (3) with 
Vsh = '^/tshi corresponding to k^ ~ 1 (see also Section 5-2b of Spitzer 1987). Recent Fokker-Planck 
calculations indicate that the rates of evaporation by two-body relaxation and gravitational shocks, 
Uf,^ and Ugh, are sometimes correlated and mutually reinforcing (Gnedin et al. 1999). The effect 
appears to be important mainly when the two rates are comparable. We neglect this complication 
and adopt Hg = 1. This and our adopted value of ^e ensure that our simple model has the cor- 
rect behavior when either two-body relaxation or gravitational shocks dominate, i.e., in the limits 
Uf,^ ^ Vgh and u^^ <^ Vg^-, and thus when the masses of clusters are either small or large (relative 
to Mp). As we show later, these limiting cases play a major role in determining the shape of the 
mass function of the clusters. 
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With the assumptions discussed above, equations (3)-(5) take the form 

M = -jXev - {l^sh + '^se)M, (6) 

where 

^ev = 269UGp)^^^mlnA, (7) 

are constants and I'se is a function of time. This has the exact solution 

M = {Mo- fiev [ exp[u,ht' + Sit')]dt'}exp[-u,ht - Sit)], (9) 

Jo 

with 

S{t) = [ Use{t')dt'. (10) 

JO 

Equation (9) can can be inverted to obtain Mo(M, t), and this can then be substituted into equation 
(2) to obtain ^{M,t) for any specified initial mass function ^q{M). 

Figure 1 shows the evolution of the mass M predicted by equations (7)-(10) for three clusters 
with different initial masses Mq on the same orbit (with the parameters specified in the caption) . In 
the early stages (t ^ 3 x 10^ yr), the mass drops approximately exponentially with time as a result 
of stellar evolution until it reaches about 60% of its initial value. Thereafter, the mass declines 
exponentially with time as a result of gravitational shocks and linearly with time as a result of 
two-body relaxation. Gravitational shocks become relatively less important with decreasing mass, 
and two-body relaxation always dominates in the late stages as the mass approaches zero. The 
evolution predicted by this simple analytical model is generally similar to that found in the more 
accurate Monte Carlo, Fokker-Planck, and N-body models, although there are differences in detail 
(see the references cited above in connection with the parameters ^e and Kg). However, even the 
most sophisticated models still involve some important idealizations, and they sometimes differ 
from each other by as much as they differ from our simple model. These comparisons indicate that 
the approximate evolution specified by equations (7)-(10) is suitable for our purposes. 

It is worth pausing here to consider separately the influence of each disruptive effect on the 
mass function. This is simplest for a set of clusters on the same galactic orbit (i.e., with the same 
values of pe^ and i^sh)- Inserting equations (9) and (10) with lyg^ = i^se = 0, pev = T^se = 0, and 
Met; = ^sh = into equation (2) gives 

ip{M, t) = ipo{M + Hevi) two — body relaxation alone, (11) 

^(M,t) = e''«''Vo(Me''«'^*) gravitational shocks alone, (12) 

'ip{M,t) = e'^(*Vo(^e'^(*)) stellar evolution alone. (13) 

For two-body relaxation alone, the masses of clusters decrease linearly with time. This flattens the 
mass function at low masses but has little effect on its shape at high masses, i.e., ^(M, t) = const 



for M < ^et)i and ip{M,t) = tpQ{M) for M > Hevt- Thus, if the mass function is initially a power 
law, ipQ^M) oc M^ with /3 < 0, it will develop a bend at M ~ ^evt- In contrast, for gravitational 
shocks or stellar evolution alone, the masses of clusters decrease exponentially or approximately 
exponentially with time. This preserves the shape of the mass function in the sense that both ip and 
M are simply rescaled by time-dependent factors. Thus, if the mass function is initially a power law, 
it will remain one at all later times: ip{M,t) oc 'tpo{M) oc M^. In the case of gravitational shocks, 
the rescaling factor increases indefinitely, whereas in the case of stellar evolution, it saturates 
at exp5'(t) ~ 1.6 for t > 3 x 10* yr. As we show later, the scaling relations derived here are 
approximately correct even when the mass function is averaged over realistic distributions of orbits 
(i.e., with different values of Hev and I'sh)- 

We now consider clusters on different orbits within a galaxy. This is assumed for simplicity to 
have a static, spherically symmetric potential ^{R), where R denotes the galactocentric radius in 
spherical (not cylindrical) polar coordinates. Thus, we can characterize each orbit by the constant 
values of the energy E and angular momentum J per unit mass and the angle 6 between the normals 
of the orbital plane and the disk. For some purposes, it is useful to reexpress E and J in terms of 
the pericenter and apocenter of the orbit, Rp and Ra, and the radius Re of a circular orbit with 
the same energy 

E = ^V;?{R,) + <^{Rc) = \{J/Rp,af + HRp,a), (14) 

where Vc is the circular velocity. We further assume that the mass distribution of the galaxy can 
be modeled as a singular isothermal sphere, with ^{R) = V^'^lni? and Vc = const. In this case, the 
mean density of a cluster is given by 



1 /3\4 
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PiE,J) = ^r- 1-ln^ ^ . (15) 



Rp 



We note that this depends mainly on the pericenter of the orbit and only weakly (logarithmically) 
on the shape of the orbit. Equation (15) is based on the formula for the tidal radius advocated 
by Innanen, Harris, & Webbink (1983), which accounts approximately for the tidal elongation of 
clusters. The mean density would be reduced by the factor (2/3)^ if the King (1962) formula for rt 
were adopted. The precise form of p is still an open issue in the case of non-circular orbits, but it is 
likely that equation (15) captures the primary dependence on E and J, which is sufficient for our 
purpose. ^ Inserting p{E, J) into equation (7) gives the rate of evaporation by two-body relaxation 
fiev as a function of E and J. 

To compute the fractional rate of evaporation by gravitational shocks Vgh as a function of 
E and J, we average the corresponding rate over R and 9, weighting by the number of clusters 
with each of these coordinates. The number of clusters at any radius is proportional to the time 



^The best observational evidence that equation (15) is at least approximately correct is the strong inverse corre- 
lation between the mean densities of globular clusters in the Milky Way and their present galactocentric distances, 
p (X ^-I'i'io.z j-ggg pjg Y q£ Innanen et al. 1983). For many position-velocity distributions, this implies a similar, if 
not identical, relation between p and Rp. 



spent there and hence inversely proportional to the radial component of the velocity there. The 
orientations of the orbits are assumed to be random. Thus, we have 

Vsh{E,J) = —j y^.j Jysh{E,J,R, 9) sin 9de. (16) 

Here, the functional dependence of I'sh on E, J, R, and 9 follows from equation (8) and the adiabatic 
correction factor A (see Appendix A and Figure 12). The radial and azimuthal periods of the orbits 
are given by 

"■f^" dR 



'" = 'k vim- <"> 

1 J f- dR ^^^^ 



Furthermore, the radial and vertical components of the velocity (relative to the disk) are given by 

Vr = [2E - {J/Rf - 2^{R)f/\ (19) 

Vz = Vt sin 9 = {J / R) sin 9, (20) 

where Vt is the transverse component of the velocity (orthogonal to the radius in the orbital plane). 
Equation (20) is valid because, just as a cluster passes through the disk, the radial and vertical 
components of its velocity are orthogonal. This simplifies our calculations substantially. Inserting 
equations (8) and (20) into equation (16), we have 

^,.iE, J) ^ if^ /" !^dR r ^<^^... (21) 

GpPrP^J^ Jr^ Vr{R) Jo sm9 

When iJ,eviE, J) and i>sh{E, J) are substituted into equations (9), (10) and (2), we obtain the mass 
function ip{M, t; E, J) of clusters on orbits specified by E and J. 

For some purposes, we are interested in how the mass function depends on position R and 
velocity V rather than energy E and angular momentum J. With this in mind, we define 
/(M, R, V, t)(iM(iR(iV as the number of clusters in the small element of mass-position-velocity 
space dMfiRdV centered on (M, R, V) at time t. We assume for simplicity that the cluster system 
is spherical and non-rotating. ^ Thus, by the Jeans theorem, we have 

/(M, R, V, t) = V(M, t- E, J)foiE, J), (22) 



■''In fact, the cluster system would develop some asphericity, even if it initially had none, as a consequence of 
the different rates at which clusters with different orbital orientations are disrupted by disk shocks. The globular 
cluster systems in many galaxies consist of both a nearly spherical, slowly rotating (halo) component and a moderately 
flattened, rapidly rotating (disk) component. In the Milky Way, about 27% of the known globular clusters are members 
of the disk component [this is the fraction of clusters in the Harris (1996, 1999) compilation with metallicities above 
the disk-halo division [Fe/H] — —0.8 specified by Zinn (1985)]. A complete analysis of the disruption of clusters would 
take these complications into account. Our simple model should provide a good approximation to the mass function 
and its dependence on galactocentric radius, the main goals of this study, because we average the rate of evaporation 
by disk shocks over orbital orientations [see equations (16)-(21)] and because, for most clusters, disk shocks are not 
the dominant disruptive process. 
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where foiE, J) is the initial distribution function, defined such that fo{E, J)dRdV is the fraction of 
clusters in the small element of position- velocity space dRdV with energies and angular momenta 
near E and J at t = 0. Equation (22) implies that the disruption of clusters does not alter their 
orbits. The mass function at the radius R = |R| is given by the integral of /(M, R, V,t) over all 
velocities; using equation (22) and evaluating the Jacobean for the transformation between (Vr, Vt) 
and {E, J), we obtain 



roo f'CO 

i;{M,R,t) = 47r / dVR / Ft/(M, R, V, t)dVT 
Jo Jo 



^ r rl^ /■^-(^'^) J^{M,t;E,J)fo{E,J) ^^ 
~ R'Mr) Jo [2E-{J/Rr-2m)V/' ' ^ ^ 

where Jm{E, R) = R[2E — 2<I>(i?)]^'^ is the maximum angular momentum of an orbit with a given 

energy E and radius R. Finally, we note that the average mass function in a volume bounded by 

radii i?i and i?2 is 

o rR2 

V;(M, t) = ^3_^3 J ij{M, R, t)R^dR. (24) 

In the following, we consider two simple models for the initial distribution function of the 
clusters. The first is the Eddington model: 

fo{E,J) cc exp(-E/a2)exp[-i(J/i?A^)']. (25) 

This has velocity dispersions (Tr = a and ax = cr[l + (-R/i?yi)^]~^ in the radial and transverse 
(i.e., orthogonal) directions, respectively, where the anisotropy radius Ra marks the transition from 
a nearly isotropic to a predominantly radial velocity distribution. In a logarithmic potential, the 
initial density profile (number of clusters per unit volume) is 

noiR) oc [1 + (R/RAfr'^R-^ (26) 

with 7 = (Vc/cr)'^. Distribution functions like the Eddington arise frequently in simulations of 
gravitational collapse, in which violent relaxation is nearly complete in the inner regions but not 
in the outer regions. The second model we consider has an initial distribution function of the form 

fo{E,J)^eM-E/a^)J-^^. (27) 

In this case, the radial and transverse velocity dispersions are gr = a and or = cr(l — /?)^'^, and 
in a logarithmic potential, the initial density profile is a power law, 

no{R) oc i?-2/3-7^ (28) 

again with 7 = (14/cr)^. We refer to this as the scale- free model. It is not clear which physical 
processes would produce a scale-free distribution function, although gravitational clustering in a 
self-similar hierarchy might do so. For our purposes, the most important difference between the 
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Eddington and scale-free models is that, in the former, the velocity anisotropy increases outward, 
whereas in the latter, it is the same at all radii. Thus, the distribution of pericenters is narrower 
in the Eddington model than it is in the scale-free model (see the formulae in Appendix B and 
Figure 13). 

Before presenting the results of our calculations, we pause here to enumerate the parameters in 
our models. These are: the escape probability parameter ^g, the circular velocity of the galaxy Vc, 
the central surface density and scale radius of the disk, T,f^{0) and Ra, the index characterizing the 
initial density profile 7, and the anisotropy radius Ra (for the Eddington model) or the anisotropy 
parameter (3 (for the scale-free model). As standard values of these parameters, we adopt ^e = 0.045, 
Vc = 220 km s~i, 2^(0) = 720 Mq pc'^, Rd = S kpc, and 7 = 2.5, Ra = 5 kpc (Eddington) or 
7 = 3.5, f3 = 0.5 (scale- free) . The first of these is the escape probability parameter in the Henon 
(1961) model, which should approximate the effective value of .^e for the pre- and post-core collapse 
evolution of clusters with a realistic spectrum of stellar masses. 

Our standard values of Vc, ^d{0), and Rj, are appropriate for the disk of the Milky Way (see 
Binney &: Merrifield 1998). For example, the standard values of Srf(O) and R^ imply that the 
surface density of the disk is 50 Mq pc~^ at the solar position, R = 8 kpc. Our standard values of 
7 were chosen so that the final density profile of the cluster system in the models would approximate 
the observed profile (see below). We have chosen the standard value of Ra to be the same as the 
median galacto centric radius of the globular clusters in the Milky Way, Rh ~ 5 kpc. This ensures 
that the velocity anisotropy at Rh in the Eddington model is the same as that at all radii in the 
scale- free model, namely aji = v2o"t- This is more radial anisotropy than appears to exist in the 
present velocity distribution of globular clusters in the Milky Way (Frenk & White 1980; Dinescu, 
Girard, & van Altena 1999). However, it is similar to the radial anisotropy of halo stars in the 
solar neighborhood (Binney & Merrifield 1998) and may be appropriate for the initial anisotropy 
of globular clusters (since clusters on elongated orbits are preferentially disrupted) . 

In the following, we explore the infiuence of different physical processes on the mass function 
by varying the parameters with respect to their standard values and comparing results from the 
Eddington and scale-free models. We calculate the mass function at times up to t = 12 Gyr, the 
age of globular clusters in the Milky Way favored in several recent studies (Gratton et al. 1997; 
Reid 1997; Chaboyer et al. 1998). 



3. Results 

The aim of this section is to present the results of our calculations and to compare them with 
observations. Figure 2 shows the mass functions of young star clusters in the merging Antennae 
galaxies and old globular clusters in the Milky Way. The former is based on deep UB VI images 
taken with the WFPC2 on HST and mass-to-light ratios derived from stellar population synthesis 
models (Zhang & Fall 1999). The latter is based on the luminosities of the 146 clusters with 
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known distances in the most recent compilation of data by Harris (1996, 1999) and M/Ly = 3. 
Both samples are believed to be reasonably complete over the ranges plotted [log(M/M0) > 3.8 
and log(M/M0) > 2.9, respectively]. If there is any incompleteness, however, it is likely that 
more clusters are missing from the low-mass ends of the functions. Here, and throughout this 
section, we have plotted the function ^(logM), the number of clusters per unit logM, against 
log(M/M0). This is related to the function ip(M) defined earlier, the number of clusters per 
unit M, by ^(logM) = (log e)''^Mip{M). The use of ^(logM) facilitates some comparisons with 
observations. The dashed parabola in Figure 2 depicts the usual lognormal mass function with 
a peak at Mp = 2 x 10^ Mq and a dispersion of cj(logM) = 0.5, corresponding to a Gaussian 
distribution of magnitudes with (My) = —7.3 and a(Mv) = 1-2 (for M/Ly = 3). 

As we have already mentioned, and as Figure 2 also demonstrates, the shapes of the mass 
functions of the young clusters in the Antennae and the old clusters in the Milky Way are very 
different. The former declines monotonically as ^(logM) oc M~^, whereas the latter rises to a 
peak at Mp ~ 2 x 10^ Mq and then declines. Another fact evident from Figure 2 is that the 
lognormal function provides a good representation of the empirical mass function (histogram) of 
globular clusters at high masses but not at low masses. For M < Mp, the observations can be fitted 
better by ^(logM) oc M, corresponding to ■ip{M) = const. This behavior in the empirical mass 
function was first shown by McLaughlin (1994). As we have pointed out here — for the first time 
we believe — the form tp{M) = const is a robust signature of evaporation by two-body relaxation. 
This behavior in the mass function can be traced to the fact that, in the late stages of evolution, 
the masses of tidally limited clusters decrease linearly with time. Thus, the time dt they spend in 
each small interval of mass dM = Mdt is independent of M, leading to ip{M) = const. ^ 

The results of our calculations are presented in Figures 3-11. In these diagrams, we plot 
the mass function at times t = 0,1.5,3,6, and 12 Gyr. The peak of ^(logM) at t = 12 Gyr 
is indicated by an upward-pointing arrow. We first explore the effects of different initial mass 
functions with the parameters fixed at their standard values. Figure 3 shows the evolution of 
the mass functions, averaged over all radii (actually, 1 < i? < 100 kpc) for the Eddington initial 
distribution function; Figure 4 shows the corresponding results for the scale-free initial distribution 
function. The four initial mass functions we consider are: (1) a power law, ipQ^M) oc M^ with 
P = —2, (2) the same power law but truncated below M = 3 x 10^ Mq, (3) a Schechter function, 
V'o(M) oc M'^exp(-M/M*) with /3 = -2 and M* = 5 x lO'^ Mq, (4) a lognormal function, 
V'o(M) oc exp{-i[log(M/Mp)/a(logM)]2} with Mp = 2 x 10^ Mq and o-(logM) = 0.5. 

In all our models, the mass function develops a peak, and at t = 12 Gyr, this is remarkably close 
to the observed peak, despite the very different initial conditions. For low-mass clusters (M < Mp), 
the disruption is dominated by two-body relaxation, which, as noted above, leads to ij){M) = const 
and ^(logM) oc M, in excellent agreement with the empirical mass function (histograms). This is 
true whether the initial mass function lies above or below the relation tp(M) = const, as illustrated 



^This is consistent with the relation dN/dtrh ~ const at smaU trh noted by Hut & Djorgovski (1992). 
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in the left- and right-hand panels of Figures 3 and 4. Even if the formation of low-mass clusters 
were suppressed entirely, they would appear later as the remnants of intermediate-mass clusters on 
their way to destruction. High-mass clusters (M > Mp) are mainly affected by stellar evolution 
and gravitational shocks. For the reasons given above, these shift \l'(logM) to smaller logM but 
leave its shape nearly invariant (when dynamical friction is neglected). Thus, the empirical mass 
function can always be matched above the peak by a suitable choice of the initial mass function, 
as illustrated in the lower panels of Figures 3 and 4, with the Schechter and lognormal functions. 
Unfortunately, neither theory nor observation provides much guidance on the form of V'o(-^) at 
high M. Because of small- number statistics, the mass function of old globular clusters in the Milky 
Way is uncertain above 10^ Mq, and that of young star clusters in the Antennae above 3 x 10^ Mq 
(see Figure 2). 

We have also computed the total number of clusters iVr(t) and their total mass Mj'(t) by 
integrating ^(M, f) and M^p{M,t) over M. Table 1 lists the present disruption rates, \Nt/Nt\ 
and \Mt/Mt\, and the present survival fractions, Nt/Nto and Mt/Mtoi for the four different 
initial conditions shown in Figure 3. In the two cases with steep initial mass functions, the survival 
fractions depend on the lower cutoff M; of the integration over M; thus, we list results for a range 
of values, M; = 1 — 10^ Mq. The other entries in Table 1 are not sensitive to Mi. We note that 
the timescales for the disruption of existing clusters, \Nt/Nt\~^ and \Mt/Mt\~^, are comparable 
to the present age (actually 1 — 4 times longer). The reason for this is that most of the clusters 
with shorter disruption timescales have already perished. The present values of INt/NtI in our 
models agree with the total disruption rate of globular clusters in the Milky Way estimated by Hut & 
Djorgovski (1992) from the empirical distribution of half-mass relaxation times {Nt = — 5ib3 Gyr~ 
and Nt = 98). Furthermore, the present values of \Mt/Mt\ are similar to the typical (median) 
disruption rates of individual globular clusters estimated by Gnedin & Ostriker (1997). We find that 
the present total number and mass of clusters, Nj- and Mt, are small fractions of their initial values, 
Nto and M^o, especially when the initial mass function rises toward low masses. The values of 
Mt/Mto, in particular, indicate that a substantial fraction of the field stars in the galactic spheroid 
could be the debris of disrupted clusters. However, even in the most extreme case considered here 
(the Schechter initial mass function with a lower cutoff at Mi = 1 Mq), the survival fraction is a 
few times larger than the ratio of the mass in globular clusters to the mass in the galactic spheroid 
(about 1%). Thus, to account for all the field stars in the spheroid by disrupted clusters, the initial 
mass function would have to rise more steeply than iPq[M) oc M^^ for M < 10^ Mq. 

In Figure 5, we plot the number density profiles of the cluster system at different times for 
the Eddington and scale-free initial distribution functions and the Schechter initial mass function. 
The profiles become flatter because clusters are destroyed faster near the galactic center, although 
a nearly steady form is reached by t ~ 1.5 Gyr. The final profiles in both models are in reasonable, 
although not perfect, agreement with the observed profile, which we have derived from the same 
compilation of data as we used for the mass function (Harris 1996, 1999). In fact, we chose the 
standard values of the parameter 7 = (Vc/a)"^, after some adjustment, to achieve this match. 
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Figures 6 and 7 show the evolution of the mass function for the same models when averaged 
over the inner and outer parts of the galaxy (i? < 5 kpc and R > 5 kpc), the boundary between 
these being close to the median radius of globular clusters in the Milky Way. In both cases, the peak 
mass is higher in the inner region. This is caused by the higher rate of two-body relaxation, resulting 
from the higher mean density, and the higher rate of gravitational shocks near the galactic center, 
the former effect being more important than the latter (see below). For the Eddington model, 
the logarithmic difference in the peak mass between inner and outer clusters is A log Mp = 0.2, 
corresponding to a difference in mean absolute magnitudes of A{Mv) = 0.5 (for constant M/Ly). 
For the scale-free model, the differences are A log Mp = 0.65 and A(My) = 1.6. For comparison, 
we find A(M\/) = 0.16 it 0.26 for the globular clusters with R < 5 kpc and R > 5 kpc in the 
Harris (1996, 1999) compilation. Thus, the radial variation of the peak mass in the Eddington 
model is consistent with observations (at the 1.3 o" level), whereas that in the scale- free model 
is not. Furthermore, the width of the mass function in the Eddington model agrees better with 
the observed one. The explanation for these differences in the mass functions can be found in 
the different radial variations of the velocity anisotropy in the two models. In the Eddington 
model, the anisotropy increases outward, causing a relatively narrow distribution of pericenters 
and hence disruption rates; in the scale-free model, the anisotropy is constant, causing a relatively 
wide distribution of pericenters and disruption rates (see Appendix B). From now on, we consider 
only models with the Eddington initial distribution function. 

The effect of changing the escape probability parameter ,^e is shown in Figure 8. As expected, 
the evolution of the mass function is slower and the peak mass is smaller for smaller ^^ and conversely 
for larger ,^e- The effective value of .^e is not known precisely, although we have argued above that 
it should be close to the Henon and our standard value (0.045) when allowance is made for both 
the pre- and post-core collapse evolution of clusters with a realistic spectrum of stellar masses. In 
fact, the similarity between the peak mass in our models and the observations indicates that the 
actual value of .^e cannot differ from our standard value by more than a factor of about two. 

In Figure 9, we show the effects of varying the velocity anisotropy radius Rj^ on the mass 
function of clusters in the inner and outer parts of the galaxy (i? < 5 kpc and R > 5 kpc). Small 
values of Ra imply predominantly radial orbits, with mostly small pericenters and hence large mean 
densities, whereas large values of Ra imply a nearly isotropic velocity distribution, with a wide range 
of pericenters and mean densities. This is why the mass function evolves faster and the peak mass 
is larger for smaller Ra- Moreover, the more radial are the orbits, the more similar are the peak 
masses in the inner and outer parts of the galaxy. For Ra = 2.5 kpc, we find A log Mp ~ 0, whereas 
for Ra = 15 kpc, we find A log Mp = 0.5. These correspond, respectively, to A(My) ~ and 1.2 
(for constant M/Ly). The former is consistent with the observed value, A{My) = 0.16 it 0.26, 
whereas the latter is probably not. Thus, a substantial degree of radial anisotropy is required in 
the initial velocity distribution for consistency with the weak radial variation in the empirical mass 
function. The present velocity distribution may in fact be nearly isotropic (Frenk & White 1980), 
but as a result of the preferential disruption of clusters on elongated orbits, the initial distribution 



15 



would have been more anisotropic. 

Figure 10 shows the effects of altering the surface density of the disk. Here, we present results 
for an exponential disk with double the standard central density, i.e., 2^(0) = 1440 MqPC^^, and 
for no disk at all. Figure 10 also indicates how our results depend on the energy-mass conversion 
factor Ks since the rate of evaporation by gravitational shocks is proportional to k^S^. As expected, 
a more massive disk causes the mass function of star clusters to evolve faster. However, in this case, 
unlike two-body relaxation, the peak mass decreases. The reason for this is that, as noted above, 
gravitational shocks shift the mass function to smaller masses while leaving its shape nearly invari- 
ant. However, gravitational shocks are less important than two-body relaxation in the disruption 
of low-mass clusters, even allowing for possible uncertainties in Hg and S^^. This conclusion is also 
supported by the N-body models of Vesperini & Heggie (1997). Thus, the peak mass and its radial 
variation and the low-mass shape of the mass function are all determined primarily by two-body 
relaxation rather than gravitational shocks, contrary to some statements in the literature. 

All our previous results were computed for a galaxy like the Milky Way, with Vc = 220 km s~^. 
It is also of interest to know how the mass function would evolve in other galaxies, with different 
circular velocities. For this purpose, we assume that the masses and sizes of galaxies scale as 
Mg (X V^ and Rg oc V^ (to satisfy the virial theorem, V^^ oc Mg/Rg). Then the mean internal 
densities of star clusters, like those of their host galaxies, scale as p oc Mg/ R? oc V^~ . Recent 
estimates of the exponent in the baryonic Tully-Fisher relation lie in the range k ~ 3-4 (Bell &: de 
Jong 2001, and references therein). For k at the lower end of the range, the mass function of star 
clusters is independent of Vc (since p is independent of Vc). For k at the upper end of the range, the 
mass function has the dependence on Vc shown in Figure 11. In this case, the peak mass decreases 
by A log Mp = 0.5 (A(My) = 1.2) as Vc increases from 110 to 440 km s~^. This is probably 
larger than allowed by observations (Harris 1991). More definitive comparisons will require better 
knowledge of the relations between M^, i?^, and Vc and possibly more complicated models for the 
internal structure of galaxies (e.g., with finite core radii). Moreover, dynamical friction, which is 
neglected in our models, may be important in galaxies with small Vc- 



4. Comparison With Previous Studies 

Several other researchers have suggested that disruptive processes would cause the mass func- 
tion of star clusters to evolve toward something like a lognormal function. Okazaki &: Tosa (1995) 
based their analysis on the survival triangle in the mass-radius (i.e., logM-logr/j) plane, defined 
by setting the characteristic timescales for evaporation by two-body relaxation and gravitational 
shocks, tev and tgh-, equal to the current time t (Fall & Rees 1977). Okazaki & Tosa assumed 
that clusters inside the triangle would exist without any loss of mass {M = Mq for tcv > t and 
tsh > t), while clusters outside the triangle would not exist at all (M = for tcv < t or tgh <t). In 
other words, the mass of each cluster was assumed to be a step function of time, with the step at 
t = niin(teD, tsh)-, rather than to decrease continuously with time as shown in Figure 1. The clusters 
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were postulated to have a power-law initial mass function and a Gaussian initial distribution of 
fj, = log(M/r)j), with f3 = 2.6 or 4.1. The mass function at later times was then obtained from this 
by integrating over rt with the limits of integration set by the survival triangle. In this approach, 
the shape of the present mass function is determined largely by the assumed shape of the initial 
distribution of fi, which is not explained. Ostriker & Gnedin (1997) followed the same approach in a 
study of the radial variation of the mass function, except that they postulated a bivariate Gaussian 
initial distribution of a; = log(M/r^) and y = logM and integrated over the survival triangle in 
these coordinates, including the side for dynamical friction. 

The evolution of the mass function of star clusters by disruptive processes was also considered 
by Elmegreen & Efremov (1997). They claimed this evolution would take the form 7/>(M, t) = 
ipQ{M)exp[—i'{M)t] with u{M) = M/M. Unfortunately, this is not correct, as one may verify 
by direct substitution into equation (1). Elmegreen & Efremov also claimed the disruption rate 
would have a strong inverse dependence on mass, u oc M^'^ with 7 ~ 2. This is based on the 
current disruption rates of surviving clusters estimated by Gnedin & Ostriker (1997). However, 
since the correlation between u and M is relatively weak, the value of 7 estimated in this way is 
quite uncertain (see Figure 2 of Elmegreen Efremov 1997). In fact, the data appear to be equally 
consistent with u oc M~^, the relation expected for the disruption of tidally limited clusters by 
two-body relaxation [see equations (6) and (7) here]. The mass function proposed by Elmegreen & 
Efremov has a peak at a few x 10^ Mq and approaches the initial mass function for higher masses, 
but its shape for lower masses, where disruption is important, differs markedly from the solutions 
presented here. 

Murali & Weinberg (1997a, b,c) studied the disruption of star clusters by two-body relaxation 
and gravitational shocks using a series of Fokker-Planck models. They followed the evolution of 
cluster systems with a halo component alone and with both halo and disk components. The initial 
distribution functions of the halo and disk components were represented by the scale-free model 
(called the Mestel sphere) and the Mestel disk, respectively, while the initial mass function was 
represented by a truncated power law. Murali & Weinberg (1997c) found that this model, with 
suitable choices of parameters, could reproduce many of the observed properties of the globular 
cluster system in the Milky Way. However, with the Fokker-Planck approach, they could only 
follow the evolution of clusters on a relatively sparse grid (5 x 4 x 5) in the variables Ra, M, and 
J/Jc- The four mass bins covered the range 1 x 10^ < M < 5 x 10^ Mq, thereby excluding the 
low-mass clusters most susceptible to disruption. In any case, since Murali & Weinberg did not 
display the mass function at later times, we cannot make useful comparisons between our results 
and theirs. 

Vesperini (1997, 1998) used analytical and semi-analytical models to study the evolution of 
the mass function of star clusters resulting from two-body relaxation, gravitational shocks, stellar 
evolution, and dynamical friction. He assumed the clusters were tidally limited and on circular 
orbits perpendicular to the galactic disk. The initial mass function was assumed to be a truncated 
power law or a lognormal function. Vesperini found in many cases that the final mass function 
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in his models resembled the empirical mass function of old globular clusters. However, for the 
truncated power law, the peak of the final mass function was well below the observed peak unless 
the truncation mass was large, i.e.. Mi > 10^ Mq. In addition, the mass functions in Vesperini's 
models have a strong dependence on galactocentric radius because, with all the clusters on circular 
orbits, no radial mixing occurs. As we have shown here, the low-mass end of the empirical mass 
function of globular clusters can be reproduced even if the initial mass function has no truncation 
(i.e., with Ml = 0). Moreover, we find that radial mixing of orbits is necessary to account for the 
weak radial dependence of the empirical mass function. Vesperini speculated that the lognormal 
mass function represented a quasi-equilibrium distribution. We find instead that the high-mass 
shape of the mass function, whatever its initial form, is approximately preserved, while the low- 
mass shape, ip{M) = const, is flatter than the lognormal function. Vesperini (2000, 2001) has 
recently used his models to predict the radial variation and dependence of the mass function of star 
clusters on the properties of their host galaxies (although the clusters are still assumed to be on 
circular orbits). 

Baumgardt (1998) studied the evolution of the mass function of star clusters resulting from 
two-body relaxation and dynamical friction but not gravitational shocks or stellar evolution. He 
assumed the clusters were tidally limited at the pericenters of their orbits and approximated their 
disruption by means of a simple analytical model. Baumgardt computed the orbits of the clusters 
numerically, with energies and angular momenta drawn from an initial distribution function similar, 
but not identical, to our scale- free model (note that 7 in his notation is 7 + 2/3 in our notation) . He 
adopted a power-law initial mass function with f3 = —2 and a broad initial distribution of half-mass 
radii. Baumgardt found that the mass function developed a peak and that, after a Hubble time, 
this coincided roughly with the peak in the empirical mass function of old globular clusters. He 
also found that the peak mass was larger at smaller radii, with AlogMp ~ 0.4 (A(My) ~ 1.0) 
for clusters inside and outside -R = 10 kpc, and he noted that this was probably incompatible 
with observations. Our results for the scale-free model generally agree with Baumgardt's, although 
detailed comparisons are difficult because his mass functions are very noisy. As we have shown 
here, the mass function in the Eddington model has a weaker radial variation than that in the 
scale-free model, in satisfactory agreement with observations. 



5. Discussion 

We have presented a series of simple, largely analytical models to compute the effects of 
disruption on the mass function of star clusters. Our models include evaporation by two-body 
relaxation and gravitational shocks and mass loss by stellar evolution. A virtue of our approach 
is that it leads to a clear understanding of how each disruptive process shapes the mass function 
of star clusters. Our goal has been to determine under what initial conditions the mass function 
evolves into a form resembling that of old globular clusters. We make two idealizations to simplify 
our calculations. First, we neglect correlations between the effects of two-body relaxation and 
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gravitational shocks. A comparison with more accurate Monte Carlo, Fokker-Planck, and N-body 
models indicates that the errors introduced by this approximation are acceptably small, especially 
for low- and high-mass clusters. Second, we assume the galactic potentials in which the clusters 
move are static and spherical. This ensures that the pericenter of each orbit remains fixed. We 
discuss below how our results would be modified if the galactic potentials were time-dependent 
and/or non-spherical. 

We find that, for a wide variety of initial conditions, the mass function in our models develops 
a turnover or peak, which, after 12 Gyr, is remarkably close to the observed peak in the mass 
function of globular clusters {Mp ~ 2 x 10^ -^o)- Below the peak, the evolution is dominated by 
two-body relaxation, and the mass function always develops a tail of the form tp{M) = const. This 
reflects the linear decrease in the masses of tidally limited clusters with time just before they are 
destroyed. The predicted form of ^(M) at and below the peak coincides well with the observed 
form. We interpret this as strong support for the idea that evaporation by two-body relaxation 
played a major role in shaping the mass function of globular clusters at low and intermediate masses 
(M < Mp). Above the peak, the evolution of the mass function is dominated by stellar evolution 
at early times and by gravitational shocks at late times (when dynamical friction is neglected). 
These processes operate at fractional rates that are independent of the masses of the clusters and 
hence tend to preserve the shape of the mass function at high masses (in a log- log plot). We also 
find that the disruption of clusters can contribute substantially to the field star population in the 
galactic spheroid if the initial mass function of the clusters rises steeply toward low masses. 

The radial variation of the mass function within a galaxy depends on the initial position- 
velocity distribution of the clusters. The reason for this is that the rate of evaporation by two- 
body relaxation and gravitational shocks depends on the orbits of the clusters, especially their 
pericenters. If most of the orbits are circular or if the velocity distribution is isotropic, the mass 
function will vary more with galactocentric radius than is observed. This variation can be reduced 
by the greater radial mixing that occurs when the velocity distribution has some radial anisotropy. 
However, to obtain a nearly uniform mass function within a galaxy, the radial anisotropy must 
increase outward, producing a distribution of pericenters that is narrower than the distribution 
of instantaneous positions of the clusters. We illustrate this by our models with Eddington and 
scale-free initial distribution functions, both of which have the same anisotropy at the median 
radius of the globular cluster system. In the former, the radial variation of the mass function is 
compatible with observations, whereas in the latter, it is not. Unfortunately, the initial position- 
velocity distribution of globular clusters is not known because most of the original clusters have 
been destroyed. However, since the destruction occurs preferentially for clusters on elongated orbits, 
the initial velocity distribution must have been more anisotropic than the present one. 

Our conclusions to this point are based on models with static, spherical galactic potentials. In 
such models, each cluster returns to the same pericenter on each revolution about the galaxy. In 
galaxies with non-spherical potentials, however, the pericenter of a cluster may change from one 
revolution to the next. This effect should help to dilute the radial variation of the mass function. 
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An even more effective mixing of pericenters and consequent homogenization of the mass function 
may occur in galaxies with time-dependent potentials. Large variations in the potential are a 
natural consequence of the formation and evolution of galaxies by hierarchical clustering. Each 
time one galaxy merges with another, the orbits of the clusters are likely to be scrambled to some 
degree by violent relaxation. In this way, the mass functions of clusters in the inner and outer 
parts of the galaxies would also be combined or averaged, erasing any prior radial variations. In 
the hierarchical picture, merging is expected to be important early in the histories of all galaxies 
and late in the histories of some galaxies. It is not clear, however, whether merging occurred 
recently enough in most galaxies to account for the observed uniformity of the mass functions of 
their globular clusters. It is doubtful, for example, that the Milky Way or Andromeda galaxies 
experienced any major mergers in the last 8 Gyr or so (otherwise, their old disks would have been 
disrupted). Nevertheless, the mass functions of globular clusters in galaxies with time-dependent 
and/or non-spherical potentials should have less radial variation than those in the idealized models 
presented here. 

Our models and some of those mentioned in the previous section have several observational 
implications. The first is that the peak of the mass function of clusters should increase with age. 
This might be observable in galaxies in which clusters formed continuously over long periods of 
time. Alternatively, the evolution of the peak mass might be observable in galaxies with bursts 
of cluster formation at different times, such as in a sequence of merger remnants. This test may 
be difficult, however, because the luminosity corresponding to the peak mass is relatively small for 
young clusters (since Mp varies more rapidly with t than M/Ly does). The second observational 
implication is that the peak mass should decrease with increasing distance from the centers of 
galaxies, unless this has been completely diluted by the mixing of pericenters discussed above. 
Searches for radial variations in the peak mass have so far been inconclusive. This test is difficult 
because the diffuse light of the galaxies also varies with radius, making it harder to find faint clusters 
in the inner regions. Finally, the strong dependence of the peak mass on the ages of clusters and 
the weak dependence on their positions within and among galaxies cast some doubt on the use 
of the peak luminosity as a standard candle for distance estimates. This method may be viable, 
however, if the samples of clusters are carefully chosen from similar locations in similar galaxies. 

Our models also have implications for attempts to understand the formation of star clusters of 
different types. The shape of the mass function above the peak is largely preserved as clusters are 
disrupted and hence should reflect processes at the time they formed. Below the peak, however, 
the shape of the mass function is determined entirely by disruption, mainly driven by two-body 
relaxation, and hence contains no information about how the clusters formed. If there were any 
feature in the initial mass function, such as a Jeans-type lower cutoff, it would have been erased. 
In our models, the only feature in the present mass function, the peak at 2 x 10^ -^o, is largely 
determined by the condition that clusters of this mass have a timescale for disruption comparable 
to the Hubble time. Thus, it is conceivable that star clusters of different types (open, populous, 
globular, etc) formed by the same physical processes with the same initial mass function and that the 
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differences in their present mass functions reflect only their different ages and local environments, 
primarily the strength of the galactic tidal field. Our results therefore support the suggestion that 
at least some of the star clusters formed in merging galaxies can be regarded as young globular 
clusters. Further investigations of these objects may shed light on the processes by which old 
globular clusters formed. 
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A. Adiabatic Correction Factors 

In this appendix, we present approximate expressions for the average adiabatic corrections to 
the energy changes within a cluster (relative to those for an impulsive response). This is important 
because the stars in the inner region of a cluster may respond nearly adiabatically to a gravita- 
tional shock, while those in the outer region may respond nearly impulsively. The local adiabatic 
correction factors for first- and second-order energy changes, Ai(x) and A2{x), are usually regarded 
as functions of the dimensionless variable x = uj{r)Tsh, where uj{r) is the orbital angular frequency 
at a radius r within the cluster, Tgh = H/Vz is the effective duration of the shock, and H is the 
scale height of the disk. Since these energy changes are caused by tidal accelerations, they are 
proportional to r^. Thus, the mass- weighted average adiabatic correction factors for first- and 
second-order changes in the total energy of a cluster are given by 



f^^ r'^Ai^2[x{r)]p{r)r'^dr 



M,2 = ^0 ;^^^-);^j2_^- - . (Ai) 



We compute Ai^2 from equation (Al) with the following approximations. Several formulae have 
been proposed for the local adiabatic correction factors (see Gnedin et al. 1999 for a summary). 
We adopt the limiting form derived by Weinberg (1994) from linear perturbation theory: 

Ai(x)=^2(x) = (l+x2)-3/2. (A2) 

This approximates the results of N-body simulations for slow shocks. Similar expressions, but with 
more negative exponents, apply for fast shocks (Gnedin & Ostriker 1999). The errors we introduce 
by using equation (A2) in both cases are acceptably small since Ai and A2 both approach 1 for 
fast shocks. 

The effects of gravitational shocks are greatest at large radii {x < 1), where the internal po- 
tential of a cluster, usually very centrally concentrated, is nearly Keplerian. Thus, we approximate 
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the orbital angular frequency by the formula 



u;(r)=(-) (^-^J . (A3) 

This is exact for circular orbits near the tidal radius. The coefficient in equation (A3) would be 
larger for radial orbits with apocenters near rt (v2 rather than 1), but the stars would then spend 
much of their time at smaller radii [rms radius = (5/8) ^'^r^]. These effects largely cancel, indicating 
that equation (A3) is approximately valid for most stars in the outer region of a cluster (where Ai 
and A2 are non- negligible) . 

Finally, in the numerical integrations of equation (Al), we approximate the density profile of 
the clusters by the simple formula 

1 / \ ^/^ 

p(r)oc- 1-- . (A4) 

In the outer region, this matches the profile of the King (1966) model, p(r) oc (1/r — l/r^)^'^ for 
r ^ rt- Furthermore, in the inner region, it has the singular behavior appropriate for core collapse 
models, p{r) oc r^^ for r ^ (Spitzer 1987). The half-mass radius of our model is r^ = O.lSrt, 
reasonably close to that of the Henon (1961) model (r/j = O.lbrt). Figure 12 shows the resulting 
average adiabatic correction factor, ^ = ^1 = A2, as a function of the dimensionless variable 

xt = io{rt)Tsh = {47TGp/3y/\H/Vz). (A5) 

This depends on E and J through p [see equations (14) and (15)] and on J, R, and 9 through Vz 
[see equation (20)]. We adopt H = 260 pc in all our calculations (Binney & Merrifield 1998). 



B. Distribution of Pericenters 

In this appendix, we derive expressions for the density of pericenter distances n{Rp) from the 
distribution function f{E, J), the density of clusters in position- velocity space with orbital energies 
and angular momenta near E and J. To do so, we introduce the auxiliary function N(E, J), defined 
such that N{E, J)dEdJ is the number of clusters with energies and angular momenta in the small 
intervals {E, E + dE) and (J, J + dJ). Then the number of clusters with pericenters in the interval 
(i?p, Rp + dRp) is given by 



/"°o / Q771 \ 

AnRln{Rp)dRp = / N[E{Rp, J), J] f — - J 



dJdRp, (Bl) 



where E{Rp, J) is the energy of an orbit with pericenter Rp and angular momentum J. Using 
equation (14) and the relation N{E, J) = 87rV/(S, J)Pr{E, J) (see equation 2-89 of Spitzer 1987), 
where Pr{E, J) is the radial period defined in equation (17), we obtain 

27]- r°° 
n{Rp) = -p3 / Jf[E{Rp, J), J]Pr[E{Rp, J), J][{J/Rpf - V,\Rp)]dJ. (B2) 

p J rip Vc 
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As in the main text, we assume the galaxy has a logarithmic potential, with Vc = const. Then, 
for the Eddington initial distribution function, we obtain 



no{Rp) oc TTV^Rp'^ / exp(-i7a;[l + {Rp/RAf])I{x){x - l)dx, (B3) 



and for the scale-free initial distribution function 
no 



/oo 
(i-x.^{-\-fx)x-'^I{x) {x-l)dx. (B4) 



Here, we have changed the variable of integration from J to x = (J/RpVc)'^ and introduced the 
parameter 7 = (Vc/cr)'^ (as before) and the function 

I{x) = / [x{z - 1) - zlnz]-^/^dz, (B5) 

where y{x) is the upper root of the equation 

a;(y- 1) = ylny. (B6) 

In Figure 13, we plot the initial density of pericenters no(-Rp), along with the initial density 
of cluster positions no(i?), for both the Eddington and scale- free initial distribution functions with 
the standard values of the parameters (7 = 2.5 and Ra = 5 kpc for the former, 7 = 3.5 and /? = 0.5 
for the latter). This shows the important result that, for the Eddington model, the distribution 
of pericenters is narrower than the distribution of cluster positions, while for the scale-free model, 
the distributions of pericenters and cluster positions have the same power-law form. Consequently, 
there is less radial variation in the disruption rates for the Eddington distribution function than 
for the scale-free distribution function. 
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Fig. 1. — Evolution of the masses of clusters predicted by equations (7)-(10) with three different 
initial masses: Mq = (1,2,4) x 10^ Mq. All three clusters have the same orbit, with Rp = A kpc, 
Ra = Q kpc, and 9 = 45°. The dotted lines show the evolution with two-body relaxation alone, the 
dashed lines with two-body relaxation and gravitational shocks, and the solid lines with two-body 
relaxation, gravitational shocks, and stellar evolution. Note that stellar evolution dominates in the 
early stages, that gravitational shocks become relatively less important as the mass decreases, and 
that two-body relaxation dominates in the late stages. 
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Fig. 2. — Histograms of the masses of young star clusters in the Antennae galaxies and old globular 
clusters in the Milky Way. The former is from Zhang &: Fall (1999); the latter is based on data 
compiled by Harris (1996, 1999). Note that the mass function of the young clusters declines 
monotonically over the entire observed range, log(M/M0) > 3.8, whereas the mass function of the 
old clusters rises to a peak and then declines. The dashed curve is the usual lognormal representation 
of the mass function with Mp = 2 x 10^ Mq and cj(logM) = 0.5, corresponding to a Gaussian 
distribution of absolute magnitudes with (My) = —7.3 and a{Mv) = 1-2 (for M/Ly = 3). Note 
that the empirical mass function of the old clusters (histogram) is shallower below the peak than 
the lognormal function. 
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Fig. 3. — Evolution of the mass function, averaged over all radii, for the Eddington initial distri- 
bution function with the standard parameters and different initial mass functions. These are: a 
power law of index f3 = —2 {top left), the same power law truncated at 3 x 10^ Mq {top right), a 
Schechter function with (3 = —2 and M^< = 5 x 10^ Mq {bottom left), a lognormal function with 
Mp = 2 X 10^ Mq and cj(logM) = 0.5 {bottom right). Each mass function is plotted at t = 0, 1.5, 
3, 6, and 12 Gyr; the arrows indicate the peak at t = 12 Gyr. The histograms depict the empirical 
mass function of globular clusters in the Milky Way. Note that the peak mass in the models is 
similar to that in the observations for the four different initial conditions. 
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Fig. 4. — Evolution of the mass function, averaged over all radii, for the scale- free initial distribution 
function with the standard parameters and different initial mass functions. These are: a power 
law of index (3 = —2 {top left), the same power law truncated at 3 x 10^ Mq {top right), a 
Schechter function with (3 = —2 and M* = 5 x 10^ Mq {bottom left), a lognormal function with 
Mp = 2 X 10^ Mq and cj(logM) = 0.5 {bottom right). Each mass function is plotted at t = 0, 1.5, 
3, 6, and 12 Gyr; the arrows indicate the peak at t = 12 Gyr. The histograms depict the empirical 
mass function of globular clusters in the Milky Way. Note that the peak mass in the models is 
similar to that in the observations for the four different initial conditions. 
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Fig. 5. — Evolution of the number density profile of the cluster system for the Eddington and 
scale-free initial distribution functions with the standard parameters and the Schechter initial mass 
function. The profiles are plotted at t = 0, 1.5, 3, 6, and 12 Gyr. To avoid a divergence in the initial 
density profile, the mass function is truncated at 100 Mq. The data points depict the empirical 
profile for globular clusters in the Milky Way. Note that the final profiles in the models are in 
reasonable agreement with the empirical profile. 
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Fig. 6. — Evolution of the mass function, averaged over inner radii (R < 5 kpc) and outer radii 
(i? > 5 kpc), for the Eddington initial distribution function with the standard parameters and the 
Schechter initial mass function. Each mass function is plotted at t = 0, 1.5, 3, 6, and 12 Gyr; the 
arrows indicate the peak at t = 12 Gyr. The histograms depict the empirical mass functions of 
globular clusters in the Milky Way in the corresponding ranges of radii. Note that the shift in the 
peak mass in the models between inner and outer radii is relatively small. 
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Fig. 7. — Evolution of the mass function, averaged over inner radii (R < 5 kpc) and outer radii 
(i? > 5 kpc), for the scale- free initial distribution function with the standard parameters and the 
Schechter initial mass function. Each mass function is plotted at t = 0, 1.5, 3, 6, and 12 Gyr; the 
arrows indicate the peak at t = 12 Gyr. The histograms depict the empirical mass functions of 
globular clusters in the Milky Way in the corresponding ranges of radii. Note that the shift in the 
peak mass in the models between inner and outer radii is relatively large. 
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Fig. 8. — Evolution of the mass function, averaged over all radii, with different values of the escape 
probability parameter ^e (as indicated), for the Eddington initial distribution function and the 
Schechter initial mass function. Each mass function is plotted at t = 0, 1.5, 3, 6, and 12 Gyr; the 
arrows indicate the peak at t = 12 Gyr. The dashed curves represent the same models with the 
standard parameters. Note that the peak mass in the models is larger for larger ^g- 
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Fig. 9. — Evolution of the mass function, averaged over inner radii (R < 5 kpc) and outer radii 
(-R > 5 kpc), with different values of the velocity anisotropy radius Ra (in kpc, as indicated), 
for the Eddington initial distribution function and the Schechter initial mass function. Each mass 
function is plotted at t = 0, 1.5, 3, 6, and 12 Gyr; the arrows indicate the peak at t = 12 Gyr. The 
dashed curves represent the same models with the standard parameters. Note that the shift in the 
peak mass in the models between inner and outer radii is larger for larger Ra- 
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Fig. 10. — Evolution of the mass function, averaged over all radii, with a more massive exponential 
disk (left panel) and with no disk (right panel), for the Eddington initial distribution function and 
the Schechter initial mass function. Each mass function is plotted at t = 0, 1.5, 3, 6, and 12 Gyr; 
the arrows indicate the peak at t = 12 Gyr. The dashed curves represent the same models with 
the standard parameters. Note that the peak mass in the models is smaller for stronger disks. 
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Fig. 11. — Evolution of the mass function, averaged over all radii, with different values of the 
circular velocity Vc (in km s~^, as indicated), for the Eddington initial distribution function and 
the Schechter initial mass function, and for /c = 4 in the relation Mg oc V^. Each mass function is 
plotted at t = 0, 1.5, 3, 6, and 12 Gyr; the arrows indicate the peak at t = 12 Gyr. The dashed 
curves represent the same models with the standard parameters. For /c < 4, the mass function has 
a weaker dependence on Vc than shown here (as explained in the text). 
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Fig. 12. — Average adiabatic correction factor for first- and second-order changes in the total energy 
of a cluster plotted against the dimensionless variable xt = uj{rt)Tsh- See Appendix A for details. 
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Fig. 13. — Initial densities of cluster positions (solid lines) and pericenters (dashed lines) for the 
Eddington and scale- free distribution functions. Note that the distribution of pericenters is narrower 
than the distribution of cluster positions for the Eddington model but not for the scale-free model. 
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Table 1. Present Total Disruption Rates and Survival Fractions 

Initial Mass Function \Nt/Nt\'' \Mt/Mt\'^ Nt/Nto)^ Mt/MtoY 

Power Law 0.064 0.022 0.00-0.02 0.08-0.19 

Truncated Power Law 0.051 0.021 0.47 0.37 

Schechter 0.079 0.049 0.00-0.01 0.04-0.10 

Lognormal 0.072 0.057 0.17 0.16 

^Rates are in Gyr~^. 

''Ranges are for a lower cutoff in iI)q{M) from Mi = 1 Mq to 10^ Mq. 



